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Abstract. This paper serves as our first effort to develop a new triangular spectral el- 
ement method (TSEM) on unstructured meshes, using the rectangle-triangle mapping 
proposed in the conference note 1241 . Here, we provide some new insights into the origi- 
nality and distinctive features of the mapping, and show that this transform only induces 
a logarithmic singularity, which allows us to devise a fast, stable and accurate numerical 
algorithm for its removal. Consequently, any triangular element can be treated as effi- 
ciently as a quadrilateral clement, which affords a great flexibility in handling complex 
computational domains. Benefited from the fact that the image of the mapping includes 
the polynomial space as a subset, we are able to obtain optimal L 2 - and ^estimates of 
approximation by the proposed basis functions on triangle. The implementation details 
and some numerical examples are provided to validate the efficiency and accuracy of the 
proposed method. All these will pave the way for developing an unstructured TSEM based 
on, e.g., the hybridizable discontinuous Galerkin formulation. 



1. Introduction 

The spectral element method (SEM), originated from Patera [28], integrates the unpar- 
alleled accuracy of spectral methods with the geometric flexibility of finite elements, and 
also enjoys a high-level parallel computer architecture. Nowadays, it has become a pervasive 
numerical technique for simulating challenging problems in complex domains [SJ [3J . While 
the classical SEM on quadrilateral/hexahedral elements (QSEM) exhibits the advantages of 
using tensorial basis functions and naturally diagonal mass matrices, the need for high-order 
methods on unstructured meshes with robust adaptivity spawns the development of triangu- 
lar/tetrahcdral spectral elements. In general, research efforts along this line fall into three 
trends: (i) nodal TSEM based on high-order polynomial interpolation on special interpo- 
lation points [5j [TTl [33] ; (ii) modal TSEM based on the Koornwinder-Dubiner polynomials 
[2T1 [TOl [19] ; and (iii) approximation by non-polynomial functions [31] [23] [4] . 

The question of how to construct "good" interpolation points for stable high-order polyno- 
mial interpolation on the triangle is still quite subtle and somehow open. The strict analogy 
of the Gauss-Lobatto integration rule on quadrilaterals/hexahedra does not exist on triangles 
[16] . though a "relaxed" rule can be constructed in the sense of [37] ■ We refer to [57J for an 
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up-to-date review and a very dedicated comparative study of various criteria for construct- 
ing workable interpolation points on the triangle. In general, such points have low degree of 
precision (i.e., exactness for integration of polynomials), and this motivates |26j the use of a 
different set of points for integration, which are mapped from the Gauss points on the reference 
square via the Duffy's transform [11] . The development of TSEM using Koornwinder-Dubiner 
polynomials as modal basis functions, generated by the rectangle-triangle mapping (i.e., the 
Duffy's transform), can be best attested to by the monograph [19] and spectral-element pack- 
age NekTar ( |http:/ /www.nektar.info/| . The analysis of this approach can be found in e.g., 
P31 HU HH [BJ- However, the main drawbacks of this approach lie in that the interpolation 
points are unfavorably clustered near one vertex of the triangle, and there is no correspond- 
ing nodal basis, making it complicate to implement. To overcome the second difficulty, a full 
tensorial rational approximation on triangles was proposed in |31j for elliptic problems, and 
extended to the Navier-Stokes problem in [3]. This approach still builds on the collapsed 
Duffy's transform with clustered grids. 

It is important to point out that the Duffy's transform not only leads to undesirable dis- 
tributions of interpolation/quadrature points, but also requires modifying the tensorial poly- 
nomial basis to meet the underlying consistency conditions (analogous to "pole conditions" 
in polar/spherical coordinates) induced by the singularity of the transform. Our mind-set is 
therefore driven by searching for a method based on a different rectangle-triangle mapping 
that can lead to favorable distributions of interpolation/quadrature points on the triangle 
without loss of accuracy and efficiency of implementation. With this in mind, we introduced 
in the conference note |24j a new mapping that pulls one side (at the middle point) of the 
triangle to two sides of the rectangle (cf. Figure [TTT1 (a)), and results in much more desirable 
distributions of the mapped LGL points (cf. Figure ITTT1 (c) vs. (d)). Moreover, this mapping 
is one-to-one. 




(a) (b) (c) (d) 

FIGURE 1.1. (a). A -H □ mapping; (b). tensorial Legendre-Gauss-Lobatto (LGL) 
points on □; (c). mapped LGL grids on A; (d). mapped LGL grids on A using the 
Duffy's transform. 



The purposes of this paper are threefold: (i) have some new insights into this mapping; 
(ii) demonstrate that the singularity of the mapping is of logarithmic type, which can be fully 
removed; and (iii) derive optimal error estimates for approximation by the associated basis 
functions. This work will pave the way for developing a new TSEM on unstructured meshes, 
which will be explored in the second part. It also brings about an important viewpoint 
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that any triangular element can be mapped to the reference square via a composite of the 
rectangle-triangle mapping and an affine mapping, and with the successful removal of the 
singularity, the triangular element can be treated as efficiently as a quadrilateral element. 
One implication is that this allows a mixture of triangular and quadrilateral elements, so one 
can handle more complex domains with more regular computational meshes, e.g., by tiling 
the triangular elements along the boundary of the domain. More importantly, for general 
unstructured triangular meshes, we can formulate the underlying variational problems using 
the recently enhanced hybridizable discontinuous Galerkin methods |5J [501 . We expect 
that the QSEM will enjoy a minimal communication between elements, and a minimum 
number of globally coupled degree of freedoms, and allow for implementing a large degree of 
nonconformity across elements (e.g., the hanging nodes and mortaring techniques). We leave 
this development to the second part after this work. 

The rest of this paper is organized as follows. In Section 2, we present some new insights of 
the rectangle-triangle mapping. In Section 3, we introduce the basis functions and the efficient 
algorithm for computing the stiffness and mass matrices with an emphasis on how to remove 
the singularity of the rectangle-triangle transform. We derive some optimal approximation 
results in Section 4, followed by numerical results on a triangle in Section 5. 



2. The rectangle-triangle mapping 

We collect in this section some properties of the rectangle-triangle mapping introduced in 
|24| . and provide some insightful perspectives on this transform. 

2.1. The rectangle-triangle mapping. Throughout the paper, we denote by 

A := {(x,y) : < x,y,x + y< l}, □ := {(£,7?) : -1 < < l}, 

the reference triangle and the reference square, respectively. The rectangle-triangle transform 
(cf. P3) T : □ -> A, takes the form 

x = i(l + 0(3-v7), y = l(3-0(l + v ), V(Oy)eD, (2.1) 

with the inversion T _1 : A — ► □ : 



U = 1 + (x - y) - VO - V? + 4(1 - x - y), 
[■q = 1 - (x -y)- yj{x - y) 2 +4(1 -x-y), 

for any (x,y) € A. It maps the vertices (—1, —1), (1, —1) and (—1, 1) of the square □ to the 
vertices (0, 0), (1, 0) and (0, 1) of the triangle A, respectively, while the middle point (1/2, 1/2) 
of the hypotenuse is the image of the vertex (1, 1) of □. In other words, this mapping deforms 
two edges (£ = 1 and r\ = 1) of □ into the hypotenuse of A, sec Figure [PI for illustration. 
Under this mapping, we have 

dx 3 — ?/ dx 1 + £ dy 1 + r\ dy 3 — £ 



<9£ 8 ' 6>77 8 ' <9£ 8 ' dn 

and the Jacobian is given by 



(2.3) 



J d(x,y) \ _ 2-t-n _ x /( x -y)2 + A(i- x - y ) _ x 



4 



M. SAMSON, H. LI & L. WANG 



For convenience of presentation, we use the handy notation: 

v = (%cy, ^ ± = (-d v ,a € ), w = (i-o^-(i-»7H, (2.5) 

where we put "~" to distinguish them from the differential operators in (x,y). Given u{x,y) 
on A, we define the transformed function: T)) = (u o T)(£, 77) = u(x,y), and likewise for 
v, etc.. Then we have 

(u,v)a = JJ u(x, y)v{x 1 y)dxdy = JJ u(£,r))v{£,r])Jd£dr). (2.6) 

Moreover, one verifies that 

V« = (d x u, d y u) = x" 1 (2(V • u) + (V T u), 2(V ■ u) - (V T u)) , (2.7) 

and 

(V«, Vu) A = ^ (V • u) (V ■ ^X-M^dr? + J ^ (VT«) (V^)x- 1 d^dr ? . (2.8) 

We observe from H2.7p - Q2.8jl that if Vu is continuous at the middle point (1/2, 1/2) of the 
hypotenuse of A, there automatically holds (note: (V T u)\n t i) = 0): 



(V- «)|(i,i) 



(du 









= 0, (2.9) 

(1,1) 

which is referred to as the consistency condition, and can be viewed as an analogy of the pole 
condition in the polar/spherical coordinates. In general, we have to build the condition (|2.9j) 
in the approximation space so as to obtain high-order accuracy, which therefore results in the 
reduction of dimension and modification of the usual basis functions (cf. [24) ) . 

One important goal of this paper is to demonstrate that this singularity can be removed, 
thanks to the observation: 

- d£dr7 = 41n2, (2.10) 

/□ 2-£-7J 

which implies that for any / e C(D), 

o 7 d t dr ) 



<4Mln2, (2.11) 



where M = maxg 77) |. In particular, the coordinate singularity can be eliminated, if / is 
a polynomial on □ (see Subsection I3.2[) . 

Now, we present other important features of this mapping. Hereafter, let / = (— 1, 1), and 
for any integer TV > 1, let Pn{I) be the set of all algebraic polynomials of degree at most TV. 
Denote by 

V N (A) :=span{xV ':0<i + j < N}, Q n (D) := (P N (I)) 2 . (2.12) 
The following property shows the correspondence between two polynomial spaces. 

Proposition 2.1. We have 

(i) ^(A)oTcSw(D). 

(ii) Qjv(D) = (Piv(A) o T) ffi x(^iv-i(A) o T). 

Here, T is the rectangle-triangle transform defined by (|2.ip . and \ = (2 — £ — 77) /2. 
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Proof. We find from {HE]) that for < i + j < N, 

i + £y/3-?7y/3-£y'/i 



x l y J = 



2 

This leads to the inclusion in (i). 
We see that for < i + j < N - 1, 

which implies x("PiV-i(A) o T) C Qat(D). 

It remains to prove Q N {D) C (V N {A) o T) © iV-i(A) o TJ, which we will show by 
induction. Firstly, by (|2.2[) . it is true for £,r), so is £77, since £ry = 5 — 4a; — 4y — 2x- Now, 
assume that it holds for {^rf with < i,j < N — 1. Then, for < i,j < N, we find 
that = ^^"V ), f r?^ = viCv 1 *' 1 ), and ^77^ = (^X^" 1 ^" 1 ) arc all of the 

form (a + bx + cy + dx)(p(x, y) + q(x, y)x)i where a, b, c, d are constants, p G Vn-i(^) and 
q G 7 , at_2(A). It is apparent that 

(a + 6x + cy + d\)(p + qx) = (a + bx + cy)p + dpx + (a + bx + cy)qx + dqx 2 

= (a + bx + cy)p + d((x — y) 2 + 4(1 — x — y))q + (dp + (a + bx + cy)q)x- 

Since (a + bx + cy)p,dx 2 q G 7-V(A) and dp, (a + bx + cy)q G 'PjV-i(A), we have 

Z N rf,ev N ,Z N V N e (^iv(A) o T) © x(Pjv-!(A) o T) , 

for all < i,j < N. This completes the induction. ■ 

In what follows, let lu > be a generic weight function on f2 = A or □. The weighted 
Sobolev space i/£ (fi) with r > is defined as in Adams [1] , and its norm and semi-norm are 
denoted by |j ■ ||r,w,n and | • |r,w,n> respectively. In particular, if r = 0, we denote the inner 
product and norm of L^(Q) by (•, and || ■ || w ,n, respectively. Moreover, if ui = 1, we drop 
it from the notation. 

Proposition 2.2. For any u G if 1 (A), we have 

^— ||V • fill 1 n + ^— IIV 1 " • Till n < IIVuL < ^-||V-u|| -in + ^IIV^-ull n , (2.13) 

2 II llx 4 11 11 x,n — 11 ha — 2 11 "x 2" Hx.n' v ' 

where x — (2 — £ — V)/^i u = uoT and the differential operators are defined in (|2.5| . 

Proof. By ([2~8]) . we have 

II l|2 11 — ||2 In — ii2 

Vli . = V ■ fi _! n + - V T fi _j n . 
II H A II llx 1 ,U 4H llx -U 

Then using the identity: 

V^fi = (1 - 00 c fi - (1 - r?)^tt = i(2 x (V ± -fi) - (£- r?)(V ■ 5)), 

|V U f A = ||V • fi|£_ 1)D + -^^(V 1 - ■ fi) - (£- »j)(V ■ fi)||^_ IiD . (2-14) 



we obtain 



As |£ — ?/| < 2, we get 

WM^ ■ fi) - (£ - »»)(V ■ «)|| J_ 1|D < 4||V X ■ G||^ n + 4||V ■ fi||"_ 1)D . 
Thus, the upper bound of (|2.13|) is a consequence of (|2.14[) . 
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It is clear that 

-4(£ - ?7)X(V ± • fi)(V • u) > -(2 X 2 |V ± • u\ 2 + 2(£ - 7,) 2 |V • u| 2 ). 

Thus, 

(2 X (V X ■ «) - (£ - r,)(V • «)) 2 > 2 X 2 |V X ■ £| 2 - (£- 77) 2 |V ■ £| 2 

>2 X 2 |V ± -u| 2 -4|V- ? i| 2 , 

which implies 

pxiV 1 - ■ u) - (£ - »y)(V • «)||*-i )P > 2|| V x • u|| 2 n - 4||V . «|| 2 _ in . 
Therefore, the lower bound of (|2.13l) follows from (|2.14l) . ■ 

Remark 2.1. We find from Proposition 12.21 that under the rectangle-triangle mapping (|2.1I) . 
the space if 1 (A) is mapped to the weighted space on □: 

:= {u G L 2 (D) : V • u G L 2 (D) , . « G £ 2 -i(n) }, (2.15) 

and vice verse. ■ 

2.2. Some new perspectives and a comparison study. Next, we have some insights of 
the rectangle-triangle mapping and compare it with the Duffy's transform |11) . 
Firstly, the transform (|2.1[) is a special case of the general mapping Tg : □ M> A : 

(..^(^ '-"-f'^' .H^), v«,„en, ,, M , 

with = 1/2. We see that this mapping pulls the hypotenuse of A into two edges of A at the 
point (9, 1 — 9). The limiting case with 9 = reduces to the Duffy's transform: To : □ n- A : 

X = ^(l+®(1-V), y=\{l + rj), V(^r/)GD, (2.17) 

with the inverse transform: T^ 1 : AkD: 

2x 

£=- 1, r) = 2y-l, V(x,y)€A. 

1 - V 

It collapses one edge, 77 = 1, of □ into the vertex (0, 1) of A. As the singular vertex corresponds 
to one edge, the Duffy's transform is not a one-to-one mapping, as opposite to (|2.1I) . This 
results in a large portion of mapped LGL points clustered near the singular vertex of A (see 
Figure O (d)). The Jacobian of ([2~TT]) - ([2~^|) is J = (1 - 77) /8, and we have 

Vu = (--—d*u, \ + ® d £ u + 2d„u) . (2.18) 

\ 1 — Tj 1 — T) J 

Different from (|2.9|1 . the corresponding consistency condition of the Duffy's transform becomes 
c>£u(£, 1) = 0. In a distinct contrast with (|2.10p . the integral 

— — d^ = oo. (2.19) 

□ 1 - V 

Consequently, the consistency condition has to be built in the approximation space, and much 
care has to taken to deal with this singularity for Duffy's transform-based methods in terms 
of implementation and analysis. 
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Secondly, the nature of the point singularity of (|2.1j) is reminiscent to that of the Gordon- 
Hall mapping |13| . which maps the reference square to the unit disc via 



V2 



y =JL^/2-e, v for,) en, 



and whose Jacobian is (2 — £ 2 — rj 2 )/ \/(2 — £ 2 )(2 — -q 2 ). It is clear that this transform induces 
singularity at four vertices of the reference square (cf. Figure 12. ip . It is worthwhile to point 
out that the collocation scheme on the unit disc using this mapping was discussed in [15| , and 
this mapping technique was further examined in [2|- 



■HMMHiMMMmmimnn 



■iiiiiiiiiiiiiiiiiiiiiiiiiiiiil 




FIGURE 2.1. Left: tensorial Legendre-Gauss-Lobatto points on the square. Right: 
the corresponding mapped LGL points on the unit disc. 



In addition, we find that the rectangle-triangle transform (|2.1[) can be derived from the 
symmetric mapping on □ : 

x = Z + V , y = £ V , V(Z,rj)en. (2.20) 

It transforms any symmetric polynomial in (£, rf) to a polynomial in (x, y), so it is referred to 
as a symmetric mapping |34J. One verifies that the image of this mapping is the curvilinear 
triangle (see Figure (b))0 

fl = { (x, y) : 1 - x + y, 1 + x + y, x 2 - Ay > 0} . 

As the symmetric mapping (|2.20[) . denoted by T : □ i— > fl, can not distinguish the images of 
(^,77) and (t?,C); it is n °t one-to-one. To amend this, one may restrict the domain of T to 
the upper triangle, denoted by A up , in □ (see Figure |2~21 (a,)), and interestingly, the square of 
maximum area contained in this subdomain is one-to-one mapped to the triangle of maximum 
area included in the curvilinear triangle Q, that is, 

f :Q : =(-1,0)X(0,1) .— > A:={(x,y):\x\<l + y<l}, (2.21) 

is a bijective mapping (see the shaded parts in Figure [2~2l (a) - (b) ) . For clarity of presentation, 
we denote the coordinate of any point in □ by (£, 77). It is clear that the reference square □ 
and □ are connected by the affine mapping: F\ : □ 1— > □, of the form (see the shaded parts 
of Figure O (a), (c)): 

C=^, <?=^. v(€,r,)en, (2.22) 

It is worthwhile to note that thanks to the symmetric mapping T : □ t— > Q, Xu 1361 discovered the first 
example of multivariate Gauss quadrature. 
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and the affine mapping: F2 : A H> A, takes the form (see the shaded parts of Figure [2721 (b). 
(d)): 

1 1 

(2.23) 



^(y + x + l), y=^{y-x + l), V(x,j)eA. 



In summary, we have □ □ 1— A i— h A . Remarkably, this composite mapping is identical 
to the rectangle-triangle mapping (|2.1[) . i.e., T = F\ o T o F 2 - 




FIGURE 2.2. (a). The reference square □, the upper triangle A up = 
{(x,y) : — 1 < x < y < 1} and the square □ (shaded), (b). The image (resp. A 
(shaded)) of the symmetric mapping T whose domain is A up (resp. □ (shaded)), 
(c). Domains obtained from □ and the upper triangle A up in (a) by the affine 
mapping Fi. (d). Domains obtained from A and f2 in (b) by the affine mapping 
F 2 . 



3. Basis functions and computation of the stiffness matrix 

We introduce in this section the modal and nodal basis functions on triangles, and present 
a fast and accurate algorithm for computing the stiffness matrix with a focus on how to deal 
with the singularity (cf. ([23u>([2~TTj) h 

3.1. Modal basis. Let / = (—1, 1) as before. We define the space 

Y N (A) = Q N (p) o T- 1 = (P N (I)) 2 o T~\ (3.1) 

which consists of the images of the tensor-product polynomials on □ under the inverse mapping 
T _1 defined in (|2.2[) . As a direct consequence of Proposition [2J] (ii), we have 

Y n (A)=V n (A)® X 'Pn-i(A), (3.2) 

where x = V ( x ~ v) 2 + 4(1 — x ~ v)i an d we recall that Pn(A) is the set of polynomials on 
A of total degree at most N. This implies that Yjy(A) contains not only polynomials, but 
also special irrational functions: for any <p G T'V-i(A). 
Define 

0o(O = ^, MO = ^-J&(Q, l<k<N-l, <MC) = i±i, (3.3) 

where J^' 1 is the Jacobi polynomial of degree k (cf. (32]). It is clear that {4>k}k=o forms a 
basis of Pn(I), and we have 

Q N (n)=spim{$ki:$ki(Z,v) = MO<f>l(v), 0<k,l<N}. (3.4) 
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It is a commonly used C°-modal basis for QSEM, which enjoys a distinct separation of the 
interior and boundary modes (including vertex and edge modes). All interior modes are zero 
on the triangle boundary. The vertex modes have a unit magnitude at one vertex and are 
zero at all other vertices, and the edge modes only have magnitude along one edge and are 
zero at all other vertices and edges. 

In view of (|3.1[) and (|3.4[) . we obtain the modal basis for Y/v(A) : 

Y N {A) = span{^, : * fcI (x,y) = ^or 1 , < k,l < N}. (3.5) 

3.2. Computation of the stiffness matrix. Though the singular integral of (|2.11l) -typc 

has a finite value, some efforts are needed to compute such integrals in a fast and stable 
manner. Next, we devise an efficient algorithm for this purpose. 

Let Lfe be the Lcgcndre polynomial of degree k, and recall that (see, e.g., [3"2"] ) 

(i - cVfcixCC) = 2ibTI^ fc - l(0 ~ Lfc+l(c) )' (3 - 6) 

(2k + l)L fc (C) = L' k+1 (Q - L' fe -i(0, (3.7) 
^(0 = ^1^(0 + ^^(0. (3.8) 



Thus, we have 



<f>'o(0 = -^MO = -4>' N (0, 0UC) = -|ifc(C). l<k<N-l. (3.9) 



(3.10) 



By ([231), (T2JJ) and (|3T4|) . 

X d x * kl = 2(4>m<t>i{ri) + MZMV)) + [(1 - - (1 - »7)0fc(O^i(»/)] ■ 

xdyVu = 2(4 (O^W + &(0$fa)) - [(i - 04(0^ fa) - (i - v)M0<t>'M] ■ 

Thanks to (|3.6j) - (|3.9|) . x^x^ki and xdy^ki can be represented by a linear combination of 
{Lk±i(^)Li±j (r])}i ,,3=0,1 • In view of this, we can evaluate the entries of the stiffness matrix by 
computing the integrals of the product of Legendre polynomials: 

s r/':=/ A V^-V^d,d, <-> g ^^f^^ ^r,:=atf. (3.11) 
Using the fact that the product L rn L n can be represented by {Lp}^^ 1 ; 

m-\-n 

= E c "" l p(0, (3-12) 

where the expansion coefficient {c™ n } can be found in e.g., [T5], we obtain 

tt ?=liE«'^ where = // d^dry. (3.13) 

Now, we describe how to compute {a P9 } in a fast and accurate manner. This essentially 
relies on the following recurrence relation. 

Lemma 3.1. We have 

Qp,q+1 ~ a p,q-l _ a p+l,g ~ a p-l,q w „ > 1 /o 

2g + 1 2p + 1 ' P ' q ~ ' 1 ' ^ 
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Proof. The statement is true for p = q > 1, since a PjP ±i = a p ±i ;P . In view of the symmetry: 



a qp , it suffices to show it holds for p > q > 1. 



We start with recalling the Lcgcndrc functions of the second kind (see, Formula (4.61.4) 
in [32]): 



(x) 



1 f 1 L n (t) 



dt, n > 1; Q (aO 



1, x + 1 
In ■ 



4 i-t ' 2 x— 1 

and the important identity (see [32] Formula (4.62.1)]): 

1 f 1 L n {x)-L n {t) 



x — t 



Vx > 1, 



dt 



(3.15) 



= ; In 



x+1 
x-l 



L n (x) - L n -i(x). 



Here, L n is the Legendre polynomial of the second kind, satisfying 



- . . 2n + 1 - , . 
Ln{x) = — xL n -i(x) 



-L n - 2 (x), n>l: L-i(x) = 0, L (x) = 1, 



(3.16) 
(3.17) 



n +1 n + 1 

which follows from (|3.16|) and [32j Formula (4.62.13)] directly. 

Using (|3.15[) - (|3.17[) and the orthogonality of the Legendre polynomials, we find that for 
p>q>l, 

1 L p (0L q ( V ) 
-i 2 - e - 77 



In ^-|) L 9 (2 - - 2L,_i(2 - 0] MO d£ 
3~0 



L,(2-0ip(S) df. 



(3.18) 



Thus, we have from (|3.7j) and integration by parts that 



2q+l 



111 



3-£nL, + i(2-0-£,-i(2-0 



1-C 



2q + l 



£ P (0d£ 



In 



3-^ i 9 +i(2 - o - £?-i(2 - o r^+i(0 - 



In 



2q + l 

3-^L, +1 (2-0-i,-i(2-0 



2p+ 1 



Working out the derivative, we obtain 



2^+1 



L p+ i(Q-£ p -i(Q 
2p+l 



p,(?+l ^p,q — l 



2q+l 



( (2-0ln(^|) 



3-£x L ?+1 (2-0-Vi(2-0 



(g + l/2)(3- 0(1-0 



^p+i(0 _ ip-i(0 



2p+ 1 



d£ 



%>+l.q — a p-l,9 


-/: 


2p+ 1 




dp+l,g dp— l 5 g 


1 


2p+ 1 


2^ 



1 L 9+1 (2 - - L,-i(2 - Vu(0 - L p-i(0 



_! (q + l/2)(3- 0(1-0 



2p+l 



dp— l,i 



j ? i i i 1 (2-oj;i i 1 (0(i-nd^ 



rl.l 



2p+l 

where we used the fact p > q and the orthogonality of Jacobi polynomials in the last step. 
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Equipped with (|3.14[) . we are able to compute {a P q}p> q accurately and rapidly. We sum- 
marize the algorithm as follows. 



Algorithm for computing {d, pq }p q=0 

1. Initialization 

(a) For p = 0, 1, ■ • ■ , 2N, compute a p o ; 

(b) For p = 1,2, ■ ■ ■ , 2N — 1, compute a p \ . 

2. For q = 2,3,--- ,N, 
For p = q, ■ ■ ■ , 2N — q, 

a>pq = a>p,q-2 + ^ , (flp+i,q-i — (3.19) 

Ap -\- 1 

Endfor of p, q. 

3. Set d pq = dq P for all < p < q < N. 



We describe below the details for computing the initial values. 
• Wc find from (|3~T8l) that 



O-pO — 



|^L p (01nj|— |d£ 



L p (01n 



3-£ 



d£ 



L p(0 ln 1 c d ^ :== a P + Pp- 

l 1 ? 



(3.20) 



It is clear that by (|3.7[) and integration by parts, 



L p (01n 



3-£ 



de = 



1 WO 

2p+lV,/_ 1 3-£ 



1 ^p-i(0 

i s-e 



It decays exponentially with respect to p, and the use of a Legendre-Gauss quadrature 
leads to an exponentially accurate approximation, since the function 1/(3 — £) is 
analytic within an ellipse (see [35]). We find from e.g., [12] that 



Pp 



Using (I3~T31) . 



Lp{i)\n — &Z 
l 1 — t 



if p = 0, 
if p > 1. 



p(p + 1) 

and the orthogonality of Legendre polynomials, we find 



'/'i 



:2-f-ij 



d£d?; — 
d£d?7 



2a p o — 



□ 2 - £ - 7? 

(g + l)«p+1.0 +P«p-1,0 



(2-QL p (Q 
□ 2-e-r? 



d£d?7 



(2 - £ - Tj)i p (0 



L P (0 dC 



d// 



p > 1. 



(3.21) 



2p+l 

We see that with an accurate computation of the initial values {a p o}i marching by 
(13.211) and (|3.19l) is expected to be stable. In Figure 13.11 wc provide a schematic 
illustration of sweeping the stencils by the algorithm. 



12 



M. SAMSON, H. LI & L. WANG 





D 


O 


o 





o 


o 




• 


• 


Diat 


;ram 



marked by "•" are marched via Steps 1-2 in the Algorithm, and those marked by 
"o" are obtained by the symmetric property in Step 3. 

Remark 3.1. We see that the rectangle-triangle mapping (|2.ip essentially induces logarith- 
mic singularity. Indeed, numerical quadrature of integrands involving a logarithmic weight 
function is of independent interest (see, e.g., [12]). ■ 



Remark 3.2. As a quick note, the mass matrix under this basis is sparse. Indeed, by (|2.6I) . 

(u, v) a = g JJ uv d£dr/ - i JJ ^uv d£d?? - JJ^WV d^dry, (3.22) 

so we claim this from (|3.8j) and the orthogonality of the Legendre polynomials. ■ 

Remark 3.3. With an additional affine mapping, any triangular element A any can be trans- 
formed to the reference square □. It is important to point out that the stiffness and mass 
matrices on A any can be precomputed in a similar fashion as above. To justify this, we con- 
sider a general triangle A any with vertices V{ = (xi,yi), i = 1,2,3. Like (|2.ip . we have the 
mapping from □ to A any : 

( \ ( , 0--00--V) , , , (1+0(3-7?) , (3-0(l + ry) 

{x,y) = [xi,yi) + {x 2 ,y 2 ) g + {X3,y 3 ) , (3.23) 

for all (£,77) G Q- A direct calculation leads to 

u )A any = -g- JJ ™ d £ d? 7 ~iqJj&v d£dr/ - — JJ rjuv d^d??, (3.24) 

and 

(Vu,Vv) Aany =A JJ (V -u)(V ■ 5)x -1 d£di7 + C JJ (V T u) (V T w)x _1 d^d7? 

- B JJ [(V"5)(V t C) + (V T u)(V-5)]x _1 d^d?j, (3.25) 

where = 2/(2 — £ — ??), the differential operators are defined in (|2.5j) . and the constants 
are given by 

F = (x 2 - an)(y 3 - 2/1) - (x 3 - x 1 )(y 2 - yi) ^ 0, 
A=((x 2 ~x 3 ) 2 + (y 2 -y 3 ) 2 )/(2F), 
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B = ((X2 - Xl ) 2 + (y 2 - yif - (x 3 - xif - (y 3 - Vif)IW, 
C = ((2H - a* - ^ 3 ) 2 + (2yi - 2/2 - 2/ 3 ) 2 )/(8F). 

In particular, if A any = A, ([XM]) and (|3~^5j) (note: 5 = 0) reduce to flU} and (|2~%]) , 
respectively. 

As with p.lOp . we find from (|3.6|) - ()3.9|) that V • and V T $fc; can be expressed in terms 
of {Lk±i(t;)Li±j(r])}ij = o t i, so the mass matrix on A any can be precomputed by the same 
algorithm described above. ■ 

3.3. Interpolation, quadrature and nodal basis. Through the general mapping (|3.23[) . 
the operations (e.g., interpolation, quadrature and numerical differentiations) on a triangular 
element can be performed on the reference square □. 

Hereafter, let {Cj}jLo be the Legendre-Gauss-Lobatto (LGL) points, i.e., the zeros of 
(1 — ( 2 )L' N (Q 7 and let {hj}jL be the associated Lagrangian basis polynomials such that 
hj G Pn(I) an d hj((k) — §kj, where Skj is the Kronecker delta. Given v G C(T), the one- 
dimensional polynomial interpolant of u is 

N 

(I C N v)(C)=J2 v ^)h j (0^PN, VCgJ. (3.26) 
3=0 

Recall that the LGL quadrature has the exactness: 

,1 n 

/ 0(CK = V0GP 2 jv_i(/), (3.27) 

where {tOj} are the LGL quadrature weights. 

Given any u G C(A), we define the interpolant of u by 

n 

(JL N u)(x,y) = (4/Mor 1 = ( (^r)te>%>«(C)fti(l)) °T~\ (3.28) 

i,j=0 

where T and T _1 are defined in (|2.1[) and (|2.2I) as before, and = r)k = Cfel^Lo- Notice 
that lL N u G Fjv(A). 

We also extend the LGL quadrature to define the discrete inner product on A as 

1 N 

(u,v) Nt A = - u{€h r )j)v{€i,'nj)x(ti,Vj)uiUj, (3.29) 

where x= (2-(-)j)/2. Asa consequence of ([23]1 . (pHT]) and (|3~T j) - (j3~2"j) . there holds 

(u,v)n,a = (u,v)a, Vu • u G y 2 jv-2(A), (3.30) 

which also holds for all it ■ u G V2N -2(A) . 

Since {h^hi}^ t _ forms the nodal basis for Qjv(D), we can obtain the nodal basis for 
Y N (A) : 

Yjv(A) = span{* feZ : $ w (x,y) = (h^oT' 1 :0<k,l<N}. (3.31) 

In view of (|3.22[) , the mass matrix under this nodal basis can be computed easily as usual by 
tensorial LGL quadrature. However, the direct evaluation of the stiffness matrix like (|3.11[) 
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is prohibitive, as there is no recursive way for the computation. In order to surmount this 
obstacle, we resort to the notion of "discrete transform" (cf. [301 )■ Like (|3.10[) . we have 

Xd^ki = 2(h' k (0hi(v) + h k (0h'M) + [(1 - OKiCMv) - (1 - ri)hk(£)h'M G Qn(P), 

xdy^u = 2(/4(£)W + hkiOhKv)) - [(i - OK(Ohi(v) - (i - v)hk(£)hfa)] e Cjv(n). 

The idea is to transform {xd x ^ki}k.i=o and {xd y ^u}k.i=o to {^(O-^j'WIij^o via a two- 
dimensional discrete transform. Then the evaluation boils down to finding {a\^ } in p. lip as 
before. 

4. Estimates of orthogonal projection and interpolation errors 

The section is devoted to error estimates of orthogonal projections and interpolations on 
triangles. These results will be essential for understanding the approximability of the basis 
functions, and provide important tools for error analysis of the TSEM for PDEs. 

4.1. Orthogonal projections. We start with considering the projection XIn : L 2 (A) — > 
Fjv(A), defined by 

(U N u-u,v) A =0, WveY N (A). (4.1) 
Theorem 4.1. For any u £ H r (A) with r > 0, we have 

\\Tl N u-u\\ A < cN~ r \u\ r ,A, (4.2) 
where c is a positive constant independent of N and u. 

Proof. We have 

J4TTI 

\\Ii N u-u\\ A = inf ||^-«a < U~u\\ A , tyeV N (A). (4.3) 

4>eY N (A) 

Now, we take t/j to be the best L 2 -approximation in "Pjv(A), denoted by tt^u. By J55J Theorem 
3-3], 

h N u-u\\ A <cN- r [ £ \\^^(ay-d x )^u\\l^^ A y /2 <cN- r \u\ r ,A, (4.4) 

ki_+k 2 +k 3 =r 

where u ' 1 > k2 < ka = x kl+k;i y k2+k3 (1 — x — y) kl+k2 is a Jacobi weight function on A. Therefore, 
the estimate (14.21) follows from (1431) and (14.41). ■ 



We now turn to the i? 1 -projection: 11^ : H 1 (A) — > Y/v(A) such that 

(V(U 1 N u-u),\/v) A + (n 1 N u-u,v) A =0, VveY N (A), (4.5) 

and the ^-projection: : H$(A) -> Y#(A) = Y N (A) n ffo( A ), defined by 

(V(U]fu-u),Vv) A = 0, VveY°(A), (4.6) 

where Hq(A) is defined as usual, i.e., the subspace of H 1 (A) with functions vanishing on the 
boundary of A. 

Theorem 4.2. For any u £ £Tq(A) n H r (A) with r > 1, we have 

||n^°«-«|| M ,A <cN^- r \u\ r , A , n = 0,l, (4.7) 

where c is a positive constant independent of u and N. It also holds for any u £ H r (A) with 
TL^u in place of 11^ it. 



A NEW TSEM I: IMPLEMENTATION & ANALYSIS 15 

Proof. Here, we only provide the proof for the projector njj°, as the estimate for LT^ can be 
obtained in a very similar fashion. 

By the Poincare inequality, we know that the semi-norm | • |i a is a norm of Hq (A). Hence, 
by the definition 



|«-HJ|. o ||i i a<c|0-u|i ) a<c||0-«||i ) a, V4>eY N (A). (4.8) 



It is known from (|3.2[) that Vn(A) c Yjv(A), so we can take <f> to be the orthogonal projection 
T^ :ffi(A)^n( 



n 1 / : ^(A) -> V% (A) = V N (A) n H£(A), defined by 



(V(7r^°«-u),Vw) A = 0, Vv£V° N (A). (4.9) 
We quote the estimate in [551 Theorem 3.4]: 

||7r^°u-u||i,A < ciV 1_r ^ £ \\d^(d y -d x )^u\\\ M)A ) 1/2 

k 1 +k 2 +k 3 =r + (4-10) 

where 

w fei,fc 2 ,fc 3 _ a ,max(fci+fe3-l,0)ymax(fc 2 +fc3-l,0) ^ ^ _ y^naax(fei+fe 2 -l,0) 

Hence, the estimate (|4.7p with /i = 1 follows from (|4.8|) and (|4.10[) . 

To show (|4.7[) with \i = 0, we use a duality argument as in [7], which we sketch below. 
Given j 6 L 2 (A) , wc consider the auxiliary problem: Find u g £ Hq(A) such that 

a(u g ,v) := (Wu g ,Wv) A = (g,v) A , V-u G Hq(A). (4.11) 



By a standard argument, we can show that this problem has a unique solution with the 
regularity ||u 9 || 2 ,a < c||g||A- 



Now, taking v = u — U]fu into (|4~lT|) . we find from (|4.6p and (|4. 7[) with /i = 1 that 



n]y°M) A | = |a(u s ,u - n^°u)| = \a{u g - XV]fu gi u- ll]fu 
< \u„ - li A r u„ L A u - II ^ u 



3 xi JV "S|i,Al JV "Ml, A 

< ciY _r |u g |2,A|w|r,A < cA^ _r ||.9i|AKi|r,A- 



Finally, we derive 



II rrbo II \(g,u-U^ u) A \ 

W-n^yU . = SUp 1|— 1| L < CN \u\ r ,A- 

o# 9 eL2(A) \\9\\a 

This completes the proof. ■ 

Remark 4.1. It is seen that benefited from the fact that 7-V(A) C Y/v(A), we are able 
to obtain the optimal error estimates directly from the available polynomial approximation 
results on triangles. ■ 

4.2. Estimation of interpolation error. Now, we estimate the error of interpolation by 
(I3.28|) on A. The estimate of the one-dimensional LGL interpolation (cf. (|3.26p ) is useful for 
our analysis (see [301 Theorem 3.44]), that is, for any v £ H r (I) with r > 1, we have 

\\I C N V - „|| < cN- r \\(l ( 2 ) {r - 1)/2 V {r) \\ LHI y (4-12) 
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Theorem 4.3. For any u 6 H r (A) with r > 2, 

\\I N u-u\\ A < cN- r B r {u), 

where 

B j u)= fM2,A + ||(^-3«) 2 «||j-i I A+||V.tt|| J -i, A) if r = 2, 

\\u[r.A + \u\r-l.A, if T > 3, 

J is the Jacobian as defined in (|2.4I) . and c is a constant independent of u and N . 



(4.13) 
(4.14) 



Proof. To this end, let I d be the identity operator. Using flU}, ([3728]) and (|4~T2|) . we obtain 
|| Hjv w - w|| A < cHi^J^u-ull,-, = c||(4 - Id)(I% - Id)u+(I}f - Id)u+(r]f - Id)u\\ D 

< c(||(4 - I d ){P N - I d )u\\ u + 11(4 - I d )u\\ a + ||(J& - I d )u\\ u ) 

< cN- 1 ^ - I d )d^u\\ a + c(||(4 - I d )u\\ D + || (J£ - /d)u|| n ) 

< d\T r (||(1 - T? 2 ) (r.-2)/2^^-l~|| n + _ e 2 )(r -l)/2^~|| n 

+ \\(l- V 2 )^)/2 d r u \\ n ). 

It remains to transform the variables (£, 77) back to (x, y) and obtain tight upper bounds of 
the right-hand side using norms of u on A. By (|2.3j) , 



c^it = — — d y u ^—{Oy -o x )u = — — d x u — [dy - d x )u, 



Thus, we have 



k=0 



1 + rj\ fl — Tj 



r-k 



(4.15) 
(4.16) 

(4.17) 



and 



|(l_^)(r-i)/a^ 



1 1 2 



i^ui 2 (i-e 2 r- i d^dTj 



where 



One verifies readily from (|2.1j) that 



<c^2JJJd r x - k (d v -d x ) k u 



k |2<2(£,r?;r,fc) 



dxdy, 



(i-er 



^(1-0(1-^=1-3:-^ 

1(1 + 0(1-7?) + 1(1 + 0(1+^ = 0:, 

l(l-0(l + r7) + l(l + 0(l+r?) = y. 
Therefore, by (|4TT5|) - (|4T2H]1 . we derive that for 2 < r < k - 1, 



(4.18) 
(4.19) 
(4.20) 



Q{£,,m r , k) 



2 k 



(1 + 0(1 - V)\ r ~ k -' /(1-€)(1-T7)X r-k-1 (1 _ ^ 
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^ „/c„ kr — k— Iz-i „ \ r — k — 1 j2 ^ „ r- 

< cx y x (1 — x — y ) J < ecu 
where we used the fact: 1 — 77 < 2 — £ — 77 = 16 J, and denoted by w 



J, 



t/(1 - x - y) 7 . 



Similarly, for 2 < r = fc, 



Q(€,V, r , k ) 



2 r 



fc «ILr-i,k,r-i,-i jA + ||(9„ -9 a ) r u" 



ro r-l,r-2,0 A 



where we used 1 — £ < 2 — £ — 77 = 16 J. Consequently, we obtain for r > 2, 
||(l-C 2 ) (r - 1)/a ^|| p 

r-l 

<c(EII 9 r fe a-^) 

< c(|u| r _i,A + Mr,A)- 

By swapping x O y and £ O 77, we get that for r > 2, 

||(1 - ^-^"lln < C(|U| P -1,A + |u| PlA ) 



(4.21) 



(4.22) 



We now turn to deal with the term II (1 — £ 2 )( r 2 ^ 2 d ri d i 



r - x u\ 



.fc=0 



r — 1\ / 1 + 77 \ 1 1 — 7] 
k 



r-fc-1 



By (I4TT61) and (l4~T71) . 



fi£"* _:l (^-^)*u 



r-l 



= E »/: r ' fc)^r fc (^ - d x ) k u + 2 W 2 (£, W, r, k)d r x - k -\d y - d x fu, (4.23) 

k=0 k=0 

where W± and W2 are polynomials of £ and 77. Thus, we have 

ll(i-a (r - 2)/2 ^r ls ir n < c E //. |g.r fc (^-^) fc H 2 (1 'j r2 d ^ 



fe=0 
r-l 



A 



+ c E //. I^r^a-^) 



^l 2(1 - a 



2Nr-2 



A 



J 



■ dxdy. 



This implies that for r > 3, 

||(l_^)(r-a)/a^-i fi || n < c (| u | r _ iiA + | u | rjA )_ 

For r = 2, we obtain from a direct calculation that 

Wd^uWa < \u\ 2 ,A + 7^W( d y ~ d x ) 2 u\\j-i^ + i||V • u||j-i,a- 
A combination of P~2Pj) - P~22"]) and (|4T23]) - (|4T25)l leads to the desired result. 



(4.24) 



(4.25) 



Remark 4.2. Like (|4.21[) . we could obtain sharper estimates with semi-norms in the upper 
bound of (|4. 1 3[) featured with the Jacobi-type weight functions w 01 ^' 1 '. 

Notice that for r = 2, the semi- norms are weighted with J -1 , as we can not factor out 
1 — £ or 1 — 77 from W\ and W2 in (|4.23p to eliminate J^ 1 . However, we point out that the 
value of JJ A J~ 1 dxdy is finite. ■ 
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5. Numerical results and concluding remarks 

In this section, we just provide some numerical results to demonstrate the high accuracy of 
the proposed algorithm for model elliptic problems on A. We also intend to compare it with 
the standard tensor-product spectral approximations on rectangles to assess the performance 
of our approach. 

Consider the elliptic equation: 

du 

-Am + 7-u = /, in A, u| ri =0, ^ = g, (5.1) 

where the constant 7 > 0, Ti is the edges x = and y = 0, T2 is the hypotenuse of A, and v 
is the unit vector outer normal to Y^. 

5.1. The scheme and its convergence. A weak formulation of (|5 . If) is to find u £ Hp (A) := 
{u £ if 1 (A) : u| ri = 0} such that 

B(« 1 »):=(Vu J V«)A+7(«.«)A = (/,»)A+7<fl > «>r a , V« G JTj^A), (5.2) 

where (-,-)r 2 is the inner product of L 2 {Y2)- It follows from a standard argument that if 
/ £ L 2 (A) and g £ L 2 ^), the problem (|5.2p admits a unique solution in Hp- (A). 

The spectral-Galerkin approximation of (|5.2|) is to find un £ Y^ 1 (A) := Yjv(A) n Hp- (A) 
such that for any vn £ Y A r 1 (A), 

B N {u N ,v N ) := (VMAr,Vujv)A + r y(uN,v N )A = (Iw/,«jv)a + (9iVn)n,t 2 , (5.3) 
where Ijv is the interpolation operator as defined in (|3.29[) , and the discrete inner product 
(g,VN)N,r2 can be defined on the quadrature rule: 

Sp 9d7= ^r[£ 1 ^ ,1)de ~/_ 1 1 ~ 7i [J2 (^(0:1) (5.4) 

where {Q, loj} are the LGL interpolation nodes and weights as before. More precisely, we 
define 

1 N 1 N 

(g,v N )N,T 2 = -j= ^2g((j, l)v N (Cj,l)cJj - -j=22g(l,Q)v N (l,Q)ujj, (5.5) 

v j=o v j=o 

where g = g a T and vn = Vn aT. 

Remark 5.1. Here, we purposely impose the Neumann boundary condition on the hypotenuse 
of A, so that the basis functions associated with this "singular" edge are involved in the 
computation. 

We reiterate that a distinctive difference with the scheme in [24] Eqn. (25)] lies in that 
the consistency condition (|2.9|) is not needed to be built in the approximation space, which 
significantly facilitates the implementation. Note that the approaches based on the Duffy's 
transform also need to modify the basis functions to meet the corresponding consistency 
condition (see, e.g., H])- ' 



To analyze the convergence of (|5.3[) . it is essential to study the approximability of the 
orthogonal projection: lY^f 1 : Hp (A) — > Y^ 1 (A), such that 

(v(n^ ri u - w ),v^) A -0, V0e^ 1(A) _ 
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Following the lines of the proof of Theorem 14.21 we find that (|4.7|) holds with U 1 ^ 1 and 
i?P (A) in place of njj and Hq(A), respectively. 

Another ingredient for the analysis is to estimate the error between the continuous and 
discrete inner products on T 2 . Using J3UJ Lemma 4.8] leads to 

\(g,v N ) N ,r 2 -(g,v N ) T2 \ <cJV-*(||(l - ^ 2 ) ( *- 1)/a S|ff(-, l)|L 3(J) ll»?J>rO, 

+ ||(l-^)(*-i)/25*5(l,.)|| L3(0 |Ml,.)|U am ). 

Then we obtain from (|4.15[) - (|4.16p and a derivation similar to the proof of Theorem 14.31 the 
following estimate: 

\(9,vn)nt 2 ~ (9,v N )r 2 \ < cN-'Wixyf-^idy - a^slUMk 

< cN-^rJvNWr,, f > 1. 

With the above preparations, we can prove the convergence of the scheme (|5.3[) by using 
Theorems I4.2fl~31 the estimate (|5.6p and a standard argument for error estimate of spectral 
approximation of elliptic problems. 

Theorem 5.1. Letu andu^ be the solutions of (|5.2p and (|5 . 3[) . respectively. Ifu G iJp^A)!! 
H r {A), f G iJ s (A) and g G # f ( r 2) wii/i r > 1, s > 2 and t > 1, tten we have 

\\u - ujv|U,a < c(i^- r |u| r ,A + N- s B s (f) + iV _t |5|t,r a ). 
where fj, = 0, 1, B s (f) is defined in (|4.14[) . and c is a positive constant independent of N and 

u,f,g- 

5.2. Numerical results. We first intend to show the typical spectral accuracy of the pro- 
posed method, so we particularly test it on (|5.1[) (with 7 = 1) with the exact solution: 

u(x,y) =e x+y - 1 sm(3xy(y-V3x/2 + V3/4)), W(x,y)eA. (5.7) 



For comparison, we also consider the standard tensor polynomial approximation of (|5.1j) on a 
square S = (0, 1/V2) 2 (note: it has the same area as A) under a similar setting, i.e., Neumann 
data on two edges x = 1/ y/2 and y = 1/ V2, and homogeneous Dirichet data on the other two 
edges. We take the exact solution: 

u(x, y) = exp( - - xj (-^= -yj^j sin (3xy(y - VZx/2 + V3/4)) , V(x,y) G 5. (5.8) 

In Figure 15.11 we plot the numerical errors of two methods, from which we observe that 
they share a very similar convergence behavior and the errors decay like 0(e~ cN ). For a fixed 
N, the accuracy of approximation on 5 seems to be slightly better than expected. We refer 
to [3J Figure 2.17] for a similar comparison of the polynomial approximations on triangles |10j 
and rectangles. Indeed, the accuracy is comparable to the existing means in [H)l |3~TI |2~4"] . 

In the second test, we choose the exact solution of (|5.1|) with finite regularity: 

u(x,y) = (l-x-y)i(e x *-l) ) V(a;,y) G A, (5.9) 
which belongs to H 3 ~ e (A) (for small e > 0). The counterpart on the square S takes the form: 

u{x,y)= (J-- x )*(-L-y)*(e?v-l), V (x,y) e S. (5.10) 
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FIGURE 5.1. Numerical errors of (|5.3|l vs. tensorial polynomial approximation 
on the square S. Left: L 2 - and L°° -errors using modal basis. Right: L 2 - and 
L°°-errors using nodal basis. 




FIGURE 5.2. Numerical errors of (|5.3[1 vs. tensorial polynomial approximation on 
the square S with solutions having finite regualrity. Left: L 2 - and L°°-errors using 
modal basis. Right: L 2 - and L°°-errors using nodal basis. 

TABLE 5.1. Comparison between the approach in [24] and the new method 



N 


Approach in [24] 


Approach in this paper 


L 2 -error 


L°°-error 


L 2 -error 


L°° error 


15 


2.866S-06 


1.018S-05 


2.349S-06 


8.281S-06 


30 


3AWE-07 


1.203S-06 


3.087S-07 


1.09LE-06 


45 


9.940S-08 


3.513S-07 


9.299S-08 


3.283S-07 



We depict in Figure [5721 the numerical errors of two approaches in log- log scale, where the 
slopes of the lines are all roughly —3 as predicted by the theoretical results (cf. Theorem 15. 1[) . 

Finally, we compare our new approach with the method in |24j (where the explicit consis- 
tency condition (|2.9|) was built in the approximation space). One can see from Table [5~T1 that 
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both approaches enjoy a similar convergence behavior. We reiterate that the new method does 
not require to modify the basis function, so with a pre-computation of the stiffness matrix, 
the triangular element can be treated as efficiently as the quadrilateral element. 

5.3. Concluding remarks. We initiated in this paper a new TSEM through presenting the 
detailed implementation and analysis on a triangle. We demonstrated that the use of the 
rectangle-triangle mapping in |24| led to much favorable grid distributions, when compared 
with the commonly-used Duffy's transform. More importantly, we showed the induced sin- 
gularity could be fully removed. It is anticipated that with this initiative, we can develop an 
efficient TSEM on unstructured meshes built on a suitable discontinuous Galcrkin formula- 
tion. This will be discussed in a forthcoming work. 

References 

[1] R.A. Adams. Sobolov Spaces. Acadmic Press, New York, 1975. 

[2] J. P. Boyd and F. Yu. Comparing seven spectral methods for interpolation and for solving the Poisson 
equation in a disk: Zernikc polynomials, Logan-Shepp ridge polynomials, Chebyshev-Fourier series, cylin- 
drical Robert functions, Bessel-Fouricr expansions, square-to-disk conformal mapping and radial basis 
functions. J. Comput. Phys., 230(4):1408-1438, 2011. 

[3] C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods. Scientific Computation. 
Springer, Berlin, 2007. Evolution to complex geometries and applications to fluid dynamics. 

[4] L. Chen, J. Shen, and C. Xu. A triangular spectral method for Stokes equations. Numer. Math.: Theory, 
Methods Appi, 4:158-179, 2011. 

[5] Q. Chen and I.M. Babuska. Approximate optimal points for polynomial interpolation of real functions in 
an interval and in a triangle. Comp. Meth. Appl. Math. Eng., 128(2):405-417, 1995. 

[6] A. Chernov. Optimal convergence estimates for the trace of the polynomial L 2 -projection operator on a 
simplex. Math. Comp., 81(278):765-787, 2011. 

[7] P.G. Ciarlet. The Finite Element Method for Elliptic Problems. North-Holland, 1978. 

[8] B. Cockburn, J. Gopalakrishnan, and R. Lazarov. Unified hybridization of discontinuous Galcrkin, mixed, 
and continuous Galcrkin methods for second order elliptic problems. SIAM J. Numer. Anal., 47(2):1319- 
1365, 2009. 

[9] M.O. Deville, P.F. Fischer, and E.H. Mund. High-Order Methods for Incompressible Fluid Flow, volume 9 
of Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 
Cambridge, 2002. 

[10] M. Dubiner. Spectral methods on triangles and other domains. J. Sci. Comput., 6(4):345— 390, 1991. 

[11] M.G. Duffy. Quadrature over a pyramid or cube of integrands with a singularity at a vertex. SIAM J. 
Numer. Anal, 19(6):1260-1262, 1982. 

[12] W. Gautschi. Gauss quadrature routines for two classes of logarithmic weight functions. Numer. Algo- 
rithms, 55(2-3):265-277, 2010. 

[13] W.J. Gordon and C.A. Hall. Construction of curvilinear co-ordinate systems and applications to mesh 
generation. Internat. J. Numer. Methods Engrg., 7:461-477, 1973. 

[14] B.Y. Guo and L. Wang. Error analysis of spectral method on a triangle. Adv. Comput. Math., 26(4):473— 
496, 2007. 

[15] W. Heinrichs. Spectral collocation schemes on the unit disc. J. Comput. Phys., 199:55-86, 2004. 

[16] B.T. Hclcnbrook. On the existence of explicit /ip-finitc clement methods using Gauss-Lobatto integration 

on the triangle. SIAM J. Numer. Anal, 47(2):1304-1318, 2009. 
[17] J.S. Hesthaven. From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex. 

SIAM J. Numer. Anal, 35(2):655-676, 1998. 
[18] E.A. Hylleraas. Linearization of products of Jacobi polynomials. Math. Scand., 10:189—200, 1962. 
[19] G.E. Karniadakis and S.J. Shcrwin. Spectral/hp element methods for computational fluid dynamics. 

Numerical Mathematics and Scientific Computation. Oxford University Press, New York, second edition, 

2005. 

[20] R.M. Kirby, S.J. Sherwin, and B. Cockburn. To CG or to HDG: a comparative study. J. Sci. Comput., 
51(1):183-212, 2012. 

[21] T. Koornwinder. Two-variable analogues of the classical orthogonal polynomials. In Theory and appli- 
cation of special functions (Proc. Advanced Sem., Math. Res. Center, Univ. Wisconsin, Madison, Wis., 



22 



M. SAMSON, H. LI & L. WANG 



1975), pages 435—495. Math. Res. Center, Univ. Wisconsin, Publ. No. 35. Academic Press, New York, 
1975. 

[22] H. Li and J. Shen. Optimal error estimates in Jacobi-weighted Sobolev spaces for polynomial approxima- 
tions on the triangle. Math. Comp., 79(271):1621-1646, 2010. 

[23] H. Li and L. Wang. A spectral method on tetrahedra using rational basis functions. Int. J. Numer. Anal. 
Model, 7(2):330-355, 2010. 

[24] Y. Li, L. Wang, H. Li, and H. Ma. A new spectral method on triangles. In Spectral and High Order 

Methods for Partial Differential Equations: Selected papers from the ICOSAHOM '09 conference, June 

22-26, Trondheim, Norway, volume 76 of Lecture Notes in Computational Sciences and Engineering, 

pages 237-246. Springer, 2011. 
[25] N.C. Nguyen, J. Peraire, and B. Cockburn. Hybridizablc discontinuous Galerkin methods. In Spectral 

and High Order Methods for Partial Differential Equations: Selected papers from the ICOSAHOM '09 

conference, June 22-26, Trondheim, Norway, volume 76 of Lecture Notes in Computational Sciences and 

Engineering, pages 63-84. Springer, 2011. 
[26] R. Pasquctti and F. Rapetti. Spectral element methods on unstructured meshes: comparisons and recent 

advances. J. Sci. Comput., 27(l-3):377-387, 2006. 
[27] R. Pasquctti and F. Rapetti. Spectral element methods on unstructured meshes: which interpolation 

points? Numer. Algorithms, 55(2-3):349-366, 2010. 
[28] A.T. Patera. A spectral element method for fluid dynamics: laminar flow in a channel expansion. J. 

Comput. Phys., 54(3):468-488, 1984. 
[29] C. Schwab, p- and hp-Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics. 

Numerical Mathematics and Scientific Computation. Oxford Science Publications, 1998. 
[30] J. Shen, T. Tang, and L. Wang. Spectral Methods: Algorithms, Analysis and Applications, volume 41 of 

Springer Series in Computational Mathematics. Springer- Verlag, Berlin, Heidelberg, 2011. 
[31] J. Shen, L. Wang, and H. Li. A triangular spectral element method using fully tensorial rational basis 

functions. SI AM J. Numer. Anal, 47(3): 1619-1650, 2009. 
[32] G. Szego. Orthogonal Polynomials, volume 23. AMS Coll. Publ., fourth edition, 1975. 
[33] M.A. Taylor, B.A. Wingate, and R.E. Vincent. An algorithm for computing Fcketc points in the triangle. 

SI AM J. Numer. Anal, 38(5):1707-1720, 2000. 
[34] H. Weber. Lehrbuch der Algebra. Erster Band, Braunschweig, 1912. 

[35] Z. Xie, L. Wang, and X. Zhao. On exponential convergence of Gcgenbaucr interpolation and spectral 
differentiation. Math. Comp., In press, 2012. 

[36] Y. Xu. Common Zeros of Polynomials in Several Variables and Higher Dimensional Quadrature. Chap- 
man & Hall / CRC, 1994. 

[37] Y. Xu. On Gauss-Lobatto integration on the triangle. SIAM J. Numer. Anal, 49(2):541-548, 2011. 



